# Required packages ####
#install.packages("quantreg")
#install.packages("doParallel")
#install.packages("snow")
#install.packages("reldist")
#install.packages("rstan")
require(quantreg)
require(doParallel)
require(reldist)

# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
MMwd<-file.path(directory, "Scripts/MM")
Datawd<-file.path(directory, "Processed")
Tabwd<-file.path(directory, "Results/Tables")

# Data ####
setwd(Datawd)
load("MM_fig.RData")

# Functions #### 
setwd(MMwd)
source("MMfun.R")

# Code begins ####
no_cores <- round(0.8*detectCores()) 
alpha<-0.05
Rep<-50000


## Change in Covariates ####
##
start_time <- Sys.time()
marginals<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=as.vector(W0$lrwage3),W1=as.vector(W1$lrwage3))
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CImarg<-apply(marginals,2,quantile,probs=c(alpha/2,1-alpha/2))

start_time <- Sys.time()
covariates<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W10e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIcova<-apply(covariates,2,quantile,probs=c(alpha/2,1-alpha/2))

start_time <- Sys.time()
coefficients<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=W10e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIcoef<-apply(coefficients,2,quantile,probs=c(alpha/2,1-alpha/2))

##
## Individual Impacts ####
### Unionization ####
start_time <- Sys.time()
unionization<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1union0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIunion<-apply(unionization,2,quantile,probs=c(alpha/2,1-alpha/2))

### Public Sector ####
start_time <- Sys.time()
pubsec<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1pubsec0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIpubsec<-apply(pubsec,2,quantile,probs=c(alpha/2,1-alpha/2))

### Manufacturing Sector ####
start_time <- Sys.time()
manuf<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1manuf0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CImanuf<-apply(manuf,2,quantile,probs=c(alpha/2,1-alpha/2))

### Race ####
start_time <- Sys.time()
race<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1race0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIrace<-apply(race,2,quantile,probs=c(alpha/2,1-alpha/2))

### Gender ####
start_time <- Sys.time()
gender<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1sex0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIgender<-apply(gender,2,quantile,probs=c(alpha/2,1-alpha/2))

### Part Time ####
start_time <- Sys.time()
part<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1part0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIpart<-apply(part,2,quantile,probs=c(alpha/2,1-alpha/2))

### Married ####
start_time <- Sys.time()
married<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1married0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CImarried<-apply(married,2,quantile,probs=c(alpha/2,1-alpha/2))

### Urban ####
start_time <- Sys.time()
urban<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1urban0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIurban<-apply(urban,2,quantile,probs=c(alpha/2,1-alpha/2))

### Experience ####
start_time <- Sys.time()
experience<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1exper0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIexperience<-apply(experience,2,quantile,probs=c(alpha/2,1-alpha/2))

### Education ####
#### High School ####
start_time <- Sys.time()
highSch<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1hsch0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIhigh<-apply(highSch,2,quantile,probs=c(alpha/2,1-alpha/2))

#### Some College ####
start_time <- Sys.time()
someColl<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1sobc0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIsomecoll<-apply(someColl,2,quantile,probs=c(alpha/2,1-alpha/2))

#### Associate Degree ####
start_time <- Sys.time()
assDeg<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1ad0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIassdeg<-apply(assDeg,2,quantile,probs=c(alpha/2,1-alpha/2))

#### College ####
start_time <- Sys.time()
college<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W1bche0e,W1=W1e)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIcollege<-apply(college,2,quantile,probs=c(alpha/2,1-alpha/2))

## Table Change in Covariates ####

Marg<-change(W0=W0,W1=W1,Boot=FALSE)
Cova<-change(W0=W10e,W1=W1e,Boot=FALSE)
Coef<-change(W0=W0e,W1=W10e,Boot=FALSE)

Union<-change(W0=W1union0e,W1=W1e,Boot=FALSE)
PubSec<-change(W0=W1pubsec0e,W1=W1e,Boot=FALSE)
Manuf<-change(W0=W1manuf0e,W1=W1e,Boot=FALSE)
Race<-change(W0=W1race0e,W1=W1e,Boot=FALSE)
Gend<-change(W0=W1sex0e,W1=W1e,Boot=FALSE)
Part<-change(W0=W1part0e,W1=W1e,Boot=FALSE)
Married<-change(W0=W1married0e,W1=W1e,Boot=FALSE)
Urban<-change(W0=W1urban0e,W1=W1e,Boot=FALSE)
Exper<-change(W0=W1exper0e,W1=W1e,Boot=FALSE)

High<-change(W0=W1hsch0e,W1=W1e,Boot=FALSE)
Soco<-change(W0=W1sobc0e,W1=W1e,Boot=FALSE)
AssDeg<-change(W0=W1ad0e,W1=W1e,Boot=FALSE)
Coll<-change(W0=W1bche0e,W1=W1e,Boot=FALSE)

prCov<-Cova/Marg
prCoe<-Coef/Marg
prGen<-Gend/Marg

prUni<-Union/Marg
prPuSec<-PubSec/Marg
prManuf<-Manuf/Marg
prRac<-Race/Marg
prGen<-Gend/Marg
prPar<-Part/Marg
prMar<-Married/Marg
prUrb<-Urban/Marg
prExp<-Exper/Marg
prHigh<-High/Marg
prSoco<-Soco/Marg
prAssDeg<-AssDeg/Marg
prColl<-Coll/Marg

A0<-actual(W0$lrwage3,Boot=FALSE)
A1<-actual(W1$lrwage3,Boot=FALSE)

T1<-t(rbind(A0,A1,Marg,Cova,Coef,
            Union,PubSec,Manuf,Race,Gend,Part,Married,Urban,Exper,
            High,Soco,AssDeg,Coll))

T2<-t(rbind(CImarg,CIcova,CIcoef,
            CIunion,CIpubsec,CImanuf,CIrace,CIgender,CIpart,CImarried,CIurban,CIexperience,
            CIhigh,CIsomecoll,CIassdeg,CIcollege))

T3<-t(rbind(prCov,prCoe,
            prUni,prPuSec,prManuf,prRac,prGen,prPar,prMar,prUrb,prExp,
            prHigh,prSoco,prAssDeg,prColl))


## Change in Coefficients ####
##
### Unionization ####
start_time <- Sys.time()
unionizationC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1union0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIunionC<-apply(unionizationC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Public Sector ####
start_time <- Sys.time()
pubseC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1pubsec0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIpubseC<-apply(pubseC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Manufacturing Sector ####
start_time <- Sys.time()
manufC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1manuf0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CImanufC<-apply(manufC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Race ####
start_time <- Sys.time()
raceC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1race0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIraceC<-apply(raceC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Gender ####
start_time <- Sys.time()
genderC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1sex0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIgenderC<-apply(genderC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Part Time ####
start_time <- Sys.time()
partC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1part0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIpartC<-apply(partC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Married ####
start_time <- Sys.time()
marriedC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1married0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CImarriedC<-apply(marriedC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Urban ####
start_time <- Sys.time()
urbanC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1urban0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIurbanC<-apply(urbanC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Experience ####
start_time <- Sys.time()
experienceC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1expe0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIexperienceC<-apply(experienceC,2,quantile,probs=c(alpha/2,1-alpha/2))

### Education ####
#### High School ####
start_time <- Sys.time()
highSchC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1high0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIhighC<-apply(highSchC,2,quantile,probs=c(alpha/2,1-alpha/2))

#### Some College ####
start_time <- Sys.time()
someCollC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1soco0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIsomecollC<-apply(someCollC,2,quantile,probs=c(alpha/2,1-alpha/2))

#### Associate Degree ####
start_time <- Sys.time()
assDegC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1ad0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIassdegC<-apply(assDegC,2,quantile,probs=c(alpha/2,1-alpha/2))

#### College ####
start_time <- Sys.time()
collegeC<-ParBootCI(R=Rep,nodes=no_cores,BootW=BootW,fun=change,W0=W0e,W1=We1coll0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
CIcollegeC<-apply(collegeC,2,quantile,probs=c(alpha/2,1-alpha/2))

## Table Change in Covariates ####

UnionC<-change(W0=W0e,W1=We1union0,Boot=FALSE)
PubSeC<-change(W0=W0e,W1=We1pubsec0,Boot=FALSE)
ManufC<-change(W0=W0e,W1=We1manuf0,Boot=FALSE)
RaceC<-change(W0=W0e,W1=We1race0,Boot=FALSE)
GendC<-change(W0=W0e,W1=We1sex0,Boot=FALSE)
PartC<-change(W0=W0e,W1=We1part0,Boot=FALSE)
MarriedC<-change(W0=W0e,W1=We1married0,Boot=FALSE)
UrbanC<-change(W0=W0e,W1=We1urban0,Boot=FALSE)
ExperC<-change(W0=W0e,W1=We1expe0,Boot=FALSE)

HighC<-change(W0=W0e,W1=We1high0,Boot=FALSE)
SocoC<-change(W0=W0e,W1=We1soco0,Boot=FALSE)
AssDegC<-change(W0=W0e,W1=We1ad0,Boot=FALSE)
CollC<-change(W0=W0e,W1=We1coll0,Boot=FALSE)


prUniC<-UnionC/Marg
prPubSeC<-PubSeC/Marg
prManufC<-ManufC/Marg
prRacC<-RaceC/Marg
prGenC<-GendC/Marg
prParC<-PartC/Marg
prMarC<-MarriedC/Marg
prUrbC<-UrbanC/Marg
prExpC<-ExperC/Marg

prHighC<-HighC/Marg
prSocoC<-SocoC/Marg
prAssDegC<-AssDegC/Marg
prCollC<-CollC/Marg

TT1<-t(rbind(A0,A1,Marg,Cova,Coef,
             UnionC,PubSeC,ManufC,RaceC,GendC,PartC,MarriedC,UrbanC,ExperC,
             HighC,SocoC,AssDegC,CollC))
             
TT2<-t(rbind(CImarg,CIcova,CIcoef,
             CIunionC,CIpubseC,CImanufC,CIraceC,CIgenderC,CIpartC,CImarriedC,CIurbanC,CIexperienceC,
             CIhighC,CIsomecollC,CIassdegC,CIcollegeC))
             
TT3<-t(rbind(prCov,prCoe,
             prUniC,prPubSeC,prManufC,prRacC,prGenC,prParC,prMarC,prUrbC,prExpC,
             prHighC,prSocoC,prAssDegC,prCollC))

Table1<-Tab(T1,T2,T3)
Table2<-Tab(TT1,TT2,TT3)

setwd(Tabwd)
rownames(Table1)[rownames(Table1)==""]<-1:length(rownames(Table1)[rownames(Table1)==""])
write.table(Table1, 'MMtabcov.csv', sep=',', row.names=T, col.names=T)

rownames(Table2)[rownames(Table2)==""]<-1:length(rownames(Table2)[rownames(Table2)==""])
write.table(Table2, 'MMtabcoef.csv', sep=',', row.names=T, col.names=T)

setwd(Datawd)
save.image("MMTables.RData")
